/*#include "JiangZao.h"
#include "fftwrap.h"*/
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Runtime.InteropServices;

namespace Droplet
{
    public class Noise
    {
        public static Warp warp = new Warp();
        public static float M_PI = (float)Math.PI;
        public struct NoiseReduction_t
        {
            public static int sample_rate;
            public static int frame_size;
            public static int ps_size;        // size of power spectrum

            public static float[] frame;
            public static float[] window;        // window function
            public static float[] ft;            // FFT coefficient
            public static float[] ps;            // power spectrum
            public static float[] Sold;          // smoothed power spectrum
            public static float[] Smin;          // minimum power spectrum
            public static float[] noise;         // noise power spectrum
            public static float[] prior;         // priori SNR
            public static float[] post;          // posteriori SNR
            public static float[] prob;          // speech/noise probability.

            public static float[] inbuf;         // overlapped add analysis
            public static float[] outbuf;
            public static float win_gain;       // gain by window function  通过汉明窗获得的增益
            public static drft_lookup[] fft_table;
            public static float[] outpcm;
            public static int adapt_count;     //适配计数器
        }
        public static void hamming(int n,float[] w)
        {
            int i;
            double tmp = 2 * M_PI / (n - 1);    /* 2*pi/(N-1) */
            float k = (float)tmp;    

            for (i = 0; i < n; i++)
            {
                tmp = 0.54 - 0.46 * Math.Cos(k * i);
                w[i] = (float)tmp;
            }
        }
        public static void RemoveDC(short[] x,int inLen)
        {
            int i,nOffset;
            float fSum = 0.0f;
            float tmp;
            short stmp;

            for(i = 0; i < inLen; i++) 
            {
                fSum += x[i];
            }
            tmp = -fSum/inLen;
            nOffset = (int)tmp;
            for(i = 0; i < inLen; i++) 
            {
                stmp = (short)nOffset;
                x[i] += stmp;

                if(x[i] > 32000) x[i] = 32000;
                else if(x[i] < -32000) x[i] = -32000;
            }
        }
        public static NoiseReduction_t[] NoiseReductionCreate(int inSampleRate,int inFrameSize)
        {
            int i,N;
            float total_win_gain;

            NoiseReduction_t[] pInst = new NoiseReduction_t[Marshal.SizeOf(400)];//(NoiseReduction_t[])calloc(sizeof(NoiseReduction_t),1);
          //  assert(0 != pInst);
            NoiseReduction_t.sample_rate = inSampleRate;
            NoiseReduction_t.frame_size = inFrameSize;
            NoiseReduction_t.ps_size = inFrameSize;

            N = NoiseReduction_t.frame_size;

            NoiseReduction_t.frame = new float[2 * N * Marshal.SizeOf(4)];   //(float[])malloc(2 * N * sizeof(float));
            NoiseReduction_t.window = new float[2 * N * Marshal.SizeOf(4)];  //(float[])malloc(2 * N * sizeof(float));
            NoiseReduction_t.ft = new float[2 * N * Marshal.SizeOf(4)];     //(float[])malloc(2 * N * sizeof(float));
            NoiseReduction_t.ps = new float[N * Marshal.SizeOf(4)];         //(float[])malloc(N * sizeof(float));
            NoiseReduction_t.Sold = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.Smin = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.noise = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.prior = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.post = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.prob = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));

            NoiseReduction_t.inbuf = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
            NoiseReduction_t.outbuf = new float[N * Marshal.SizeOf(4)];//(float[])malloc(N * sizeof(float));
          /*  memset(pInst->inbuf,0,N*sizeof(float));
            memset(pInst->outbuf,0,N*sizeof(float));
            memset(pInst->prob,0,N*sizeof(float));*/

            NoiseReduction_t.fft_table = Warp.spx_fft_init(2 * N);
            NoiseReduction_t.outpcm = new float[2 * N * Marshal.SizeOf(4)];//(float[])malloc(2 * N * sizeof(float));

            hamming(2 * N, NoiseReduction_t.window);
            total_win_gain = 0;
            for(i = 0; i < 2*N; i++)
                total_win_gain += NoiseReduction_t.window[i];
            NoiseReduction_t.win_gain = N / total_win_gain;
            NoiseReduction_t.adapt_count = 0;      //适配计数器初始
            return pInst;
        }
        public static void NoiseReductionDestroy(NoiseReduction_t[] inHandle)
        {
            if (inHandle == null)
                return;
        
            if (NoiseReduction_t.frame != null)
                NoiseReduction_t.frame = null;
            if (NoiseReduction_t.window != null)
                NoiseReduction_t.window = null;
            if (NoiseReduction_t.ft != null)
                NoiseReduction_t.ft = null;
            if (NoiseReduction_t.ps != null)
                NoiseReduction_t.ps = null;
            if (NoiseReduction_t.Sold != null)
                NoiseReduction_t.Sold = null;
            if (NoiseReduction_t.Smin != null)
                NoiseReduction_t.Smin = null;
            if (NoiseReduction_t.noise != null)
                NoiseReduction_t.noise = null;
            if (NoiseReduction_t.prior != null)
                NoiseReduction_t.prior = null;
            if (NoiseReduction_t.post != null)
                NoiseReduction_t.post = null;
            if (NoiseReduction_t.prob != null)
                NoiseReduction_t.prob = null;
            if (NoiseReduction_t.inbuf != null)
                NoiseReduction_t.inbuf = null;
            if (NoiseReduction_t.outbuf != null)
                NoiseReduction_t.outbuf = null;
            if (NoiseReduction_t.fft_table != null)
                NoiseReduction_t.fft_table = null;
            if (NoiseReduction_t.outpcm != null)
                NoiseReduction_t.outpcm = null;
            inHandle = null;
        }
        public static void FrequencyAnalysis(NoiseReduction_t[] inHandle,short[] x)       //X是原始数据，inHandle是main中创建的噪音缩减对象
        {
            int i;
            int N = NoiseReduction_t.ps_size;
            float[] ps = NoiseReduction_t.ps;
            float[] ft = NoiseReduction_t.ft;

            RemoveDC(x,N);   //删除直流

            for(i = 0; i < N; i++)
                NoiseReduction_t.frame[i] = NoiseReduction_t.inbuf[i];
            for(i = 0; i < N; i++)
                NoiseReduction_t.frame[i + N] = x[i];
            for(i = 0; i < N; i++)
                NoiseReduction_t.inbuf[i] = x[i];

            for(i = 0; i < 2*N; i++)
                NoiseReduction_t.frame[i] *= NoiseReduction_t.window[i];

            Warp.spx_fft(NoiseReduction_t.fft_table, NoiseReduction_t.frame, ft);     //快速傅里叶变换

            ps[0] = ft[0]*ft[0];
            for(i = 1; i < N; i++)
                ps[i] = ft[2*i-1]*ft[2*i-1] + ft[2*i]*ft[2*i];
        }
        public static void NoiseEstimation(NoiseReduction_t[] inHandle)
        {
            float[] Sold = NoiseReduction_t.Sold;
            float[] Smin = NoiseReduction_t.Smin;
            float[] noise = NoiseReduction_t.noise;
            float[] ps = NoiseReduction_t.ps;
            int N = NoiseReduction_t.ps_size;
            float beta    = 0.96f;    // look-ahead factor,control adaptive time   预取系数，控制自适应时间
            float gamma   = 0.998f;
            int i;

            //若第一次适配
            if (1 == NoiseReduction_t.adapt_count)     
            {
                float[] post = NoiseReduction_t.post;  //之后的
                float[] prior = NoiseReduction_t.prior;  //之前的

                // assume the first frame is noise     假设第一个frame是噪音
                for(i = 0; i < N; i++)
                {
                    Sold[i]  = ps[i];
                    Smin[i]  = ps[i];
                    noise[i] = ps[i];
                    post[i]  = 1.0F;
                    prior[i] = 0;
                }
            }

            for(i = 0; i < N; i++)
            {
                float Sp;
                // smooth power spectrum
                double tmp = 0.8f*Sold[i] + 0.2*ps[i];
                Sp = (float)tmp;
                // update noise   更新噪音数组
                if(Sold[i] <= Sp)
                    noise[i] = gamma*noise[i] + (1.0f-gamma)/(1.0f-beta)*(Sp-beta*Sold[i]);
                else
                    noise[i] = Sp;

                Sold[i] = Sp;
            }
        }
        public static void AudioEnhancement(NoiseReduction_t[] inHandle)
        {
            int i;
            float[] ft = NoiseReduction_t.ft;
            float[] Sold = NoiseReduction_t.Sold;
            float[] noise = NoiseReduction_t.noise;
            float[] prior = NoiseReduction_t.prior;
            float[] post = NoiseReduction_t.post;
            int N = NoiseReduction_t.ps_size;

            ft[0] = 0;    // ignore DC component  忽视直流

            for(i = 1; i < N; i++)
            {
                float gain;   //增益
                float offset;    //偏移量
                float priorSNR,postSNR;     //之前、之后信噪比
                float alpha = 0;
                // compute priori SNR and post SNR
                postSNR = Sold[i]/noise[i];
                priorSNR = postSNR-1;
                offset = (priorSNR-prior[i])/postSNR;
                alpha = 1/(1+offset*offset);
                float max = priorSNR > 0 ? priorSNR : 0;
                prior[i] = alpha*(4.0f/M_PI)*prior[i]+(1-alpha)*max;

                // compute weiner gain
                gain = prior[i] / (prior[i] + 1);

                prior[i] = priorSNR;
                post[i] = postSNR;

                // apply gain
                ft[2*i-1] *= gain;
                ft[2*i]   *= gain;
            }
        }
        //滤除噪音后数据又会放会入ioPcm，write到新文件中
        public static void NoiseReductionProcess(NoiseReduction_t[] inHandle,short[] ioPcm)
        {
            int i;
            int N = NoiseReduction_t.ps_size;

            //数据处理
            NoiseReduction_t.adapt_count++;

            FrequencyAnalysis(inHandle,ioPcm);    //频率分析
            //Console.WriteLine("频率分析成功！！！！");
            NoiseEstimation(inHandle);     //噪音估算
            //Console.WriteLine("噪音估算成功！！！！");
            AudioEnhancement(inHandle);       //语音改善
            //Console.WriteLine("语音改善成功！！！！");

            //处理完成恢复数据格式
            Warp.spx_ifft(NoiseReduction_t.fft_table, NoiseReduction_t.ft, NoiseReduction_t.outpcm);    //快速傅里叶逆变换
            // overlap and add
            for (i = 0; i < NoiseReduction_t.frame_size; i++)
                ioPcm[i] = (short)(NoiseReduction_t.outbuf[i] + NoiseReduction_t.outpcm[i]);    //j将处理后的数据放入ioPcm
            for (i = 0; i < NoiseReduction_t.frame_size; i++)
                NoiseReduction_t.outbuf[i] = NoiseReduction_t.outpcm[NoiseReduction_t.frame_size + i] * NoiseReduction_t.win_gain;

            // limit value of adapt_count     适配计数器极限值
            if (NoiseReduction_t.adapt_count > 16000) NoiseReduction_t.adapt_count = 2;
        }
    }
}


/*
 * Hamming
 *                        2*pi*k
 * w(k) = 0.54 - 0.46*cos(------), where 0 <= k < N
 *                         N-1
 *
 * n window length
 * w buffer for the window parameters
 */














